dblogr/

GxE with Mixed Models

A tutorial on how to run asses genotype by environment interactions in R


Introduction

Visual

# devtools::install_github("derekmichaelwright/agData")
library(agData)
library(lme4)      # For mixed-effects models
# Prep data
myCaption <- "derekmichaelwright.github.io/dblogr/ | Data: AGILE"

Data

# Prep data
dd <- read.csv("data_mixed_models.csv") %>%
  mutate(REP = DTM - DTF,
         VEG = DTF - DTE)
DT::datatable(dd)

GxE

GxE and Plotting Functions

myGxE <- function(xx, myTraits = c("DTE","DTF","DTS","DTM","VEG","REP")) {
  # Initialize an empty list to store variance component data
  var_list <- list()
  blup_list <- list() #store blups values
  # Loop over each trait
  for(i in myTraits) {
    # Define the model formula with ENV as a fixed factor
    formula <- as.formula(paste(i, "~ Expt + (1 | Name) + (1 | Name:Expt) + (1 | Expt:Rep)"))
    # Fit the mixed-effects model
    model <- lmer(formula, data = xx)  # Replace with your actual data
    # Extract variance components for random effects
    var_i <- VarCorr(model) %>% as.data.frame(vc)
    # Calculate total phenotypic variance 
    var_total <- sum(var_i$vcov) 
    # Calculate proportions as percentages
    var_i$Proportion <- (var_i$vcov / var_total) * 100
    #
    var_i <- var_i %>% mutate(Trait = i) %>%
      select(Component=grp, Trait, Proportion, vcov, sdcor)
    # Store results in the list
    var_list[[i]] <- var_i
  }
  # Combine all trait data into one data frame
  var_table <- do.call(rbind, var_list) %>%
    mutate(Trait = factor(Trait, levels = myTraits),
           Component = factor(Component, 
              levels = c("Residual", "Expt:Rep", "Name:Expt", "Name"),
              labels = c("Residual + Env", "Rep:Env", "G x E", "Genetic")))
  # Output
  var_table
}
gg_GxE <- function(xx) {
  myColors <- c("darkred", "purple4", "steelblue", "darkgreen")
  # Create the stacked bar plot
  mp <- ggplot(xx, aes(x = Trait, y = Proportion, fill = Component)) +
    geom_col(position = "stack", color = "black", alpha = 0.7) +
    scale_fill_manual(values = myColors) +
    theme_agData_col() +
    labs(title = "Proportion of Phenotypic Variance by Trait",
         y = "Proportion of Total Variance (%)",
         caption = myCaption)
}

All Environments

# Save variance table 
var_table <- myGxE(dd)
write.csv(var_table, "gxe_All.csv", row.names = F)
# Plot
mp <- gg_GxE(var_table)
ggsave("gxe_All.png", mp, width = 6, height = 4)

Temperate

# Save variance table 
var_table <- myGxE(dd %>% filter(MacroEnv == "Temperate"))
write.csv(var_table, "gxe_Temperate.csv", row.names = F)
# Plot
mp <- gg_GxE(var_table)
ggsave("gxe_Temperate.png", mp, width = 6, height = 4)

Mediterranean

# Save variance table 
var_table <- myGxE(dd %>% filter(MacroEnv == "Mediterranean"))
write.csv(var_table, "gxe_Mediterranean.csv", row.names = F)
# Plot
mp <- gg_GxE(var_table)
ggsave("gxe_Mediterranean.png", mp, width = 6, height = 4)

South Asia

# Save variance table 
var_table <- myGxE(dd %>% filter(MacroEnv == "South Asia"))
write.csv(var_table, "gxe_SouthAsia.csv", row.names = F)
# Plot
mp <- gg_GxE(var_table)
ggsave("gxe_SouthAsia.png", mp, width = 6, height = 4)


dblogr/


© Derek Michael Wright